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BASIC NOTATION 

f - the gravitational constant, 

E - the mass of the planet, 

M - the mass of the satellite, 
t - dynamical time, 

X, y, z - planetocentric coordinates of the satellite, 
u = X + y 

s = X - y tT^, 

P - the planetocentric distance of the satellite, 
w = 1/p, 

P(j - the variational part of p, 

Wjj - the variational part of w, 

n - the mean daily sidereal motion of the satellite, 

a - the mean semi-major axis of the satellite defined by means of the Kepler relation 

= f (E + M), 

a - the mean semi-major axis defined as the constant factor attached to the variational 
solution, 

e - the constant of the eccentricity of the satellite, 

7 - the sine of one half of the orbital inclination of the satellite relative to the orbit 

of the Sun, 

c(n - n' ) - the anomalistic frequency of the satellite, 

- the part of c independent from e, e< , y and a , 
g(n - n' ) - the draconitic frequency of the satellite, 

o; - the part of g independent from e, e <, y and a , 

J = exp (n - n') t 
D = Cd/dC, 

e ' - the eccentricity of the solar planetocentric orbit. 
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a' - the semi-major axis of the solar orbit, 

n ' - the mean daily motion of the sun in its orbit around the planet, 

m = n'/ (n - n') , 

a = a/a' - the parallactic factor, 

n - the disturbing function. 
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ON A TRANSFORMATION OF THE DIFFERENTIAL EQUATIONS 

OF THE LUNAR THEORY 


by 

Peter Musen 

Goddard Space Flight Center 


INTRODUCTION 

In recent times electronic machines have been extensively applied to develop analytical theories 
of satellites. The works by Barton (1966) and Broucke (1968) on the analytical expansion of the lunar 
disturbing function, by Broucke and Garthwaite (1969), Charpont and Mangency (1965), Danby et al 
(1965) on the manipulation of the literal series on the machine, by Deprit (1969) and Kamel (1969) on 
the application of Lie series to the elimination of the short period terms from a Hamiltonian, a beau- 
tiful work by Deprit and Henrard on the duplication of the Delaunay theory (1970), and also the work 
by Elmabsout (1970) on a semi-numerical method in the lunar theory deserve special attention. 

In this article, we suggest a transformation of the differential equations of the Hill (1878, 1886) 
and Brown (1896-1899) theory of satellites to a form integrable by the process of iteration. We re- 
tain the classical form of expansion of rectangular coordinates into trigonometric series in four argu- 
ments and powers of the satellite's mean orbital eccentricity e, sine of one half of the mean orbital 
inclination 7 , eccentricity of the solar orbit e', and parallactic factor a = a/a'. The planetocentric 
orbit of the sun is assumed to be elliptic. We differ from the classical approach in the form of the 
differential equations and method of integration. 

Hill and Brown use the method of undetermined coefficients to integrate the differential equa- 
tion for terms with a given characteristic. They form a set of equations satisfied by the coefficients 
of the expansion of coordinates, and solve these equations by iteration. This way of determining the 
coefficients is especially convenient if the computation is performed by using a desk calculator. 

At present, the explicit formulation of these equations in detail cannot be considered necessary. 

We now have analytic programming languages which permit the process of iteration to be ap- 
plied directly to the differential equations. To make this application possible, we transform Hill's 
differential equations for the complex coordinates to a form containing on the left-hand side a linear 
differential operator with constant coefficients, and on the right-hand side the small parameter m^ 
as a factor. Such a form allows the application of the method of iteration and guarantees its fast 
convergence with the accepted numerical accuracy. 

Besides the differential equations for the coordinates, we also establish (Andoyer, 1926) the 
differential equation for the parallax. This equation is redundant, but its use provides a good check 
of our results. In particular, we found that by the combined use of the differential equations for the 
coordinates and for parallax, we can determine the periodic variational solutions by the method of 
iteration, thus by-passing Hill's equations for the coefficients of the variational terms. The same 
set of equations can be used to determine the terms containing the first power of the eccentricity 
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as a factor, providing we know the main part Cp of the anomalistic frequency c. We can derive c^ 
by supression of certain terms in the expansion of the parallax. 

We suggest the use of the expansion of Andoyer (1926) for the main parts of the anomalistic 
and draconitic frequencies. Using the iterative procedure for the terms of higher order, we deter- 
mine, at each iteration step, the coefficients of terms with a given characteristic from the set of 
pairs of linear equations with two unknowns. These equations are automatically formed and solved 
by the machine. 

We derive the higher order terms in the expansion of c and g in powers of e^, e' ^ and 
by applying the averaging process to the equations of motion written in a suitable form. A close 
look at the integration procedure of the Hill-Brown theory shows that integration consists of applying 
the inverses of some qua^atic differential operators with constant coefficients to the right-hand sidef 
of the differential equations of motion. These operators represent the first terms of the expansions 
in powers of e^, e' ^ ^ 2 and of the original differential operators standing in the left -hand -sides 
of the differential equations of motion. Thus, expanding the differential and integral operators in 
powers of squares of the basic parameters constitutes an essential part of the theory. These 
expansions are also important in forming the equations satisfied by terms of higher orders in the 
development of c and g. Applying the expansion of the double differential operator to the terms 
proportional to e or , , and suppressing the secular terms in the parallax and the z- coordinate, can 
yield the equations for terms of the expansion of c and g. 

The process of integration in the Hill- Brown theory is connected with the literal expansion of 
the coordinates, frequencies and operators in powers of basic parameters. However, the literal 
expansion in m in the Hill-Brown theory is not as essential as inthe theories based on removal of 
short period terms by means of canonical transformations. Thus, to simplify the computations, 
we replace the problem of analytical convergence by the problem of numerical convergence and 
suggest, together with Brown, substituting the numerical value of m from the outset for any particular 
satellite. 

In solving the differential equations we must derive the expansions of w^,(the cube of the 
parallax) and , (the square of the distance of the satellite from the planet) in powers and products 
of the perturbations of parallax w = 1/p. This is accomplished by expanding the corresponding 
Taylor operator in powers and products of the perturbations in w and into series in so-called 
Faa-de-Bruno differential operators (Faa-de-Bruno, 1855). These last operators are associated 
with the formation of higher derivatives of a function of a function. The Taylor operator thus expandet 
is then applied to and p ^ . The differential, integral and Taylor operators in their expanded 
form, together with the extraction operator, constitute a complete set of operators to be used in 
deriving a literal theory in rectangular coordinates by iteration. 

The extraction operator performs the extraction of the terms of the desired type from the 
right-hand sides of the differential equations of motion. It is imbedded into analytical programming 
languages. In the application to the satellites of outer planets, series will, normally, never be very 
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long. The number of numerically significant terms with a given characteristic will be of order few 
tens and the number of iterations will not be much larger than five. 

THE BASIC DIFFERENTIAL EQUATIONS 

The starting point of our exposition is Hill's system of differential equations of motion of 
the satellite: 



D2u + 2mDu + — m2(u + s)-ftrw3u = - _ 

2 3 s 

(1) 


D2 s - 2m Ds + — m2 (u + s') - << s = - -5J1, 

2 B u 

(2) 

where 

(D2 _ m 2 ) z _ K w3 z = , 

2 3 z 

(3) 


f (E 4 M) ,, ,2 3 

X ~ - (1 4- a 

(n - n')2 

(4) 


The form of expansion of the coordinates is 

P = 1 j i 2 Pj, p' = 1 j ' i 4- 2 p', q = I k I * 2 Qj, r r I 2 i 1 + 2 Tj, 

where p^ , p' , Qj , are positive integers. The monomial \ = eP e'p’ a' is the characteristic 
of the term in the expansion; p + p' + q + r is its order. The linear factor a is defined in such a 
manner that the constant term in the variational part of u a. is equal to 1. 

We shall reduce the differential equations (1) - (3) to a form suitable for the computation 
of terms with the prescribed characteristic using the iterative procedure. This form is convenient 
for programming and computation on electronic machines. We introduce dimensionless coordinates, 
replacing u,s,z byua,sa,az, respectively. Equations (1) - (3) preserve their form after this 
change in the meaning of notations. The definition of « , however, will not be as in (4), but will 
become 

« = (1 + . (5) 


Substituting the variational solution 

Ug = Ul + m2 P (?2^]_ 

Sg = [1 - m2 p 
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where P is the Laurent series in ^2 ^ into (1) - (2), we deduce for the new value of k defined by 

( 5 ), 


K = Kg + S Af, 


where 


1 o ^2 

K„ = l+ 2m+ — m"', 

0 2 

and S K is of the fourth order relative to m . We obtain the numerical value of S /< by iteration, during 
the process of deriving the variational solution . The condition now must be imposed that 

the constant terms in and are equal to 1. 

We introduce a set of partial differential operators, 

A, A', r, and Y, 


acting on the expression of the type 


^2i + jc + kg + 


m 


in the following manner: 


A0=j 4>, 4> = } ' <P> r <p = k (p, X <p = 2 i (p- 


Evidently 


D = Y + cA + gr + m A'. 


The expansions 


= ^0 + ^^2 + + • • • ’ g = go + §2 + S4 + • • • > 

where and ^ 2 j are the polynomials of the jth order in squares of the basic parameters e, e ', 
y, and a, are not initially known and must be determined step by step. Consequently, neither the 
operator D nor its inverse D'* can be used. The theory of Hill and Brown is formulated so that 
not D and D"', but the operator 


Do=T + 


Cg A + gg r + m A ' 


and its inverse are used throughout. We have the expansion 

D = Dg + + D 4 + . . . . 
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where 


^^2. = ■ 82s ■"- 

s = 1, 2, 3 


We have 

Dg 0 = (2 i + Cp j + gg k + ra j '■) 0, 

°2s <^= (j k 4 >, 

Dgi(jb=(2i + Cg j +ggk +m j 

Following Andoyer (1926) we set u = S (1 + p), s = (1 + q). Equations (1) and (2) then 

become 


[D^ + 2 ( 1 + m) D + /<g] p + »<g 


(1 +p) (1 +q) = 


a n 

3 q 


( 6 ) 


[D^ - 2 (1 + m) D + K-g] q 


x -/<w3(l+q)+-|m2 
“ 2 dp 


(7) 


2 m + ^ 

2 


We rewrite them in a form more convenient for applying the process of iteration, i,e., 

D^+2(l .m)D.|.g] p,|.gq = P-||. ( 8 ) 

D2 -2(1 +m-)D+|.gJ q+|/<gP = Q-|^. (9) 


where we set 


P = - 2 + <1^ + '<0 


(v/3 _ 1) (1 + p-) (p + q) 


+ (1 + p) S 


(10) 


Q = --| m2 (2 ^ ^ 


(W3 _ I) (1 


+ ql (p + q)J 


+ w^ (1 +q) S K. 


( 11 ) 
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The quantities p and q are complex conjugates. Consequently, only one of the equations (10) and (11) 
can be considered. We can set n = 0 if we are interested only in terms depending upon e and 7 . 

Taking 

- (1 + p) (1 ^ q) 4- z 2 (12) 

into account, we can put (10) into the form 

p = _im2 (1 + q) (13) 

* -0 -f- ) - (w^ - 1) p - I (p q T z 2 ) 

convenient for using the iterative procedure. We now split p, q, w and P into the variational part, 
and the part which contains all terms with positive characteristics, and set 


P = So ~ 'S’ (I = Vq V> W = Wg ^ y (14) 

P - (Hq + S k) + H. (15) 

The zero- subscript refers to the variational part. From (14) we deduce 

p2 = p2_2y,A, (16) 

+ 3 Wp y ^ B, (11) 


where A and B are of order not smaller than two. When the terms with the given characteristic 
are to be determined, then the terms with the same characteristic as contained in A and B must 
always represent the "cross-effects" between the terms of some lower characteristics, namely, 
between the terms which were computed at a previous step. 

The analytical representation of these cross-effects will be discussed later in this exposition. 
From (12), (14) and (16) we deduce the relation 

A = 2 Y + (1 + T7g) f + (1 + 77 - p + z2). (18) 
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Substituting (14) into (13) we obtain 


c 


0 “ 



(1 + no> 


(1^0 * «'o - I) + 


-f^o ^o] 


[(1 


or, decomposing into factors, 


Ho r^d 




- (™0 - !'> ^0 -f^o "'oj 

+ [w3 (1 . - 1] 8 


This representation of is suitable for the iterative procedure, because the three terms on the 
right-hand side are of second, fourth, and sixth order in m , respectively. The differential equations 
for Cg and are 


D2 . 2 (1 * D, * I 


3 

■ 2 '"o "'o 


b K , 


(19) 


where 




- 2 (1 . m) D„ 


^0 


3 

■ 2 “"o 


( 20 ) 


Ho = H„. 
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II 


We actually need only one of these equations. They are a paraphrase of the differential equations 
[D2 + 2 (1 + m) Dq + Kq] _ k w 3 (1 + +|m2 5-2 (1 + 7^^) = 0, (21) 

[P2 _ 2 (1 + m) Dg + /<o] + /<o - « w3 (1 + Tjol + |-m2 J2 ^ j. ^ q, (22) 


easily obtainable from the general form (6) and (7), but (19) and (20) are in a form more convenient 
for using an analytical programming language. 

In a similar manner we deduce from (13), taking (12), (16) and (17) into account, that 

= = ('^ '*'0 - '^0 - I '^0 %) ^ 



^0 


T 3 [ K w5 (1 _ Kg] X' 

[3w5^y.B(l +^g+S')] [A - (f T, * z2)]. 


Eliminating y from the last equation by means of (18), we obtain H in the form suitable for actual 
computations, 

+ + (23) 


where we set 


M D 



K w3) 


= - — [k (v»2 _ 1) j. S k] , 
2 ° 


(24) 
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^ - f -'^0 “ '*'0 


S= I « w= (1 + g-p) [A - (f 77 + z^-)] 


+ K [3 w5 + B (1 + g„ +01- 

Taking into account (8), (13), (15) and (23) - (26), we conclude that the aggregate of all terms with 
positive characteristics satisfies the differential equation 

[d2 + 2 (1 + m') D + I Kg] f + I T7 = M ^ + N 77 + (s - 

We can rewrite the last equation in the form 


D 2 + 2 (1 Dg 




- d fl 

T)=Mg+Nr7+S - - — 
3 77 


(27) 


- [2 (Og i 1 + m-) S D + (5 D)2] 6 


where 


Sn = D-Dg = A5c + rsg, 

S C = C - Cq. g = g - gg- 

This form is suitable for computations and for deriving the terms with a given characteristic (or 
order), because the coefficients M and n contain m ^ as a factor. 

In the process of deriving the terms with the prescribed characteristic \ , the corresponding 
parts contributed by s and [2(Dg + 1 + m) § D + (sD)^] ^ represent the cross-action of terms 
with characteristics lower than A . The expansion of s starts with the second-order terms. The 
expression 2(Dg + 1 + m) S D ^g contains the terms c^ , c^, . . . of the expansion of c linearly and 
is closely associated with the formation of the linear equations for their determination. The parts 
and 77 ^ , which are proportional to the first power of e, satisfy the differential equation 
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(28) 


+ 2 a + m) D„ +1 K J + 1 7,^ = M + N V^, 

[d^ - 2 a + in') Dg 4 I /< J ■*• I '"O M ■^e + ^e- (28') 

To derive the expansion of the parallax one can use the equation established by Andoyer (1926): 


K W + 



3 m ^ 


(29) 


irti2 (1 
2 


m D”*) u • 


3 2/1 

+ _m2 (1 


m D ') - C + 4> = 0, 


where 


<i> = n + 


If sn 

— u 

2\ 3 u 


3 n 

3 s ' 


B n \ 

3 z / 


■!• m D 


B n B a B n\ 

3u Bs Bi-C'j 


C , the constant of integration, is developable into a series 

= ^0 ^2 + *-4 ■ ■ ■ ’ 

where C^. are homogeneous polynomials of jth degree in , e' ^ and . 

Substituting (14) into (29) we can split the differential equation into the equation for the 
variational part in parallax 

K Wg + +^°o) Po 

+ | n >2 (1 4 tnDg ') (1 + ? g )"] 


4 |.m2 (1 - m DgM (1 + = ^0’ 


and the differential equation for the aggregate of all terms with positive characteristics 
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(D^ _ K + 2 m2) Y = K (w3 _ 1) y 


( 29 ") 


+ (m2+iD2) 


A — 3 m ^ z ^ 


+ 3 m2 (1 + m D->) [c ^2 (i + + i 

+ 3 m2 (1 _ m D"‘) [ ^"2 (1 + 77^^ + -jTj) 77] - S C + >I>. 

We do not follow in the footsteps of Andoyer and do not use the upper equation for the determination 
of y . The expansion of y can be derived from (18) directly, without additional integration. We 
suggest, however, the use of the differential equation for the effects proportional to e, 

(D 2 _ k + 2 m2) y^ = « (w2 - 1 ) (30) 

+ 3m2 (1 xmD-l) (1 + ^„) 

+ 3 m2 (1 _ m (1 + vj . 

either to derive the value of , if one prefers not to use an analytic equivalent to Hill's determinan- 
tal equation, or to check the expansions if necessary. 

We can expand y into a series 


y = Vi + V2 + Xj + • • ■ . 

where are homogeneous polynomials of kth degree, in e, e', y and a, with the coefficients which 
are power series in C" , and i"' . We also need the expansion of and into a series in 
powers and products of y^ , y^ , y^ , . . .. 

Designating by T the Taylor operator 


T = expf w® y ^ 

we have 

w2 = T w2, 
P^^Tpl 

The operator T can be expanded into the series 
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(31) 


in’-..-'*.-*.. ■■ 


k = 1 , 2 , • • . j. 


j = 1, 2, 3, 


where I is the identity operator and T. are Faa-de-Bruno (1855) differential operators. They 
are the homogeneous polynomials of kth degree relative to the operators 


and satisfy the recursive relations 


or, in matrix form, the relation 


s = 1, 2, . . . 


ST. ^ 

a J “ cr, k “ 1 


Sj 0 0 0 0 


Sj 0 0 0 


T 22 0 0 0 


Sj Sj 0 0 0 


T,, 0 0 


T T 0 0 

32 33 " 


S 3 Sj 0 0 


• ^3 T 3 , T 33 0 


= T T T 0 

42 43 '44 


S4 83 S^ S, 0 


^4 "^42 ^43 ”^44 


It follows from the homogeneity of T, that 




^ = '"f T. ^ (yj. Yj. ■ ■ . Vu'» 


thus the symbols T ^ shall designate the polynomials T. ^ (vi > V 2 > • • • > Xk compute them 

step by step we can use either (32) or (32') by replacing S^ by , By applying the general formula 


f (w) = T f (Wjj) 
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■Ill Mini II III mil II II null i i i 


■I II I 


II II II 


to the particular case 


f (w) - p'^, 


and taking (16) and (17) into consideration, we obtain 

B . . B 3 . . 


where 


A. = + 3 T. ^ - 4 T. 


and 


0 S.2 - '♦'"S ‘j,3 + • • ■ + 




(- iv (j + n T. 


= 3 wj T. 3 + w® T. 3. for j ^3. 


(33) 

(34) 
(34') 


The actual expansion of the polynomials Tj.u in powers of ^ , C' , and S'" , and powers of the 
basic parameters, should be performed on the electronic machine . However, for general orienta- 
tion and for the purpose of checking, we give here the explicit form of these polynomials up to the 
eighth order: 


’•'3.2 = 2 Yi Y2- 
^3,3= 

3'4,2 ^ 2yj V 3 + 

3 ^ 4, 3 “ 3 yi Y2’ 

3'4,4 = yi’ 

'^ 5 , 2 = 2 Yi Y4 + 2 Vj Y 3 - 

^ 5.3 = 3 Yi Y 3 + 3 Xj 

3 ^ 5 ,4 “ ^ Yi Y2’ 

■ 35 , 5 = vf. 

'36.2 = 2 Y, Y5 + 2 Yj Y4 + y|. 

3'6. 3 = 3 + 6 Yj Xj Y 3 + y3, 

' 36.4 = ^ V? Y 3 + 6 Y? yl< 

3'e.5 = 5 yJ Y2> 

*6.6 - Yl' 

13 
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’^7 . 2 ^ 2 Vj Yg ^ 2 yj + ^ ^3 ><'4’ 

’’'7 . 3 " ^ Vs + ® 5^1 >^2 ><'4 + 3 ><'1 ^3 + 3 ^2 Vs’ 

■’^7,4 = X'l V 4 " Y 2 V3 + >'1 4 ’ 

T7.5-5yty3-i°v?v^ 

■^ 7.6 = ® X"! >^ 2 ’ 

^7.7^ Vl. 

'^ 8,2 ^ ^ ^1 ^7 ^ ^ Y 2 ^6 ^ Ys ^S Y 4 ' 

'^8,3 " ^ Vf Xe ^ Xi V 2 Yj + 6 Yj Y 3 Y^ + 3 + 3 Y 2 y^ 

"^8,4 = Yi Vs -f 12 Yi Y 2 ^4 + ^ y? Ya + 1^ Yi y| Y3 + Y 2 

Tg.s = 5 y1 y4 ^ 20 y/ y 2 v 3 + 10 y? y|- 

"^s.e = 6 Y^ yj + 15 y\ yj, 

1'8,7 = ^ Y? yj’ 

T 8 , 8 = y 1 - 

VARIATIONAL SOLUTION 

We derive Hill's variational solution of the lunar theory by applying the method of iteration 
directly to the system 


[CDg + 1 + (/<g + m2)j + 1 k ^ = Hq + S « , 


( 19 ') 


- R„ 


where 


3 2 

:„ = m'^ 

0 2 


(1 + V 


( 29 ’) 


+ '<o{ K - 1'>^ [3 pi h + y + K 

+ (W^ - 11 $0 - I ?o '?o} 

+ [w^ (1 - IJ i K, 
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Ro = 




+ |m2 (1 +mD-»-) a ^ ^^^->2] 
-|m2 (1 -mD->) [r^* (1 4 77^)2] . 


The values of , -q^, and > ^ from the previous iteration can be substituted into and . 

We start the iteration with some suitable approximation to ^ and . In the absence of 

better information, we can even begin with the values 

" 4 = 0, Cp r,^ 0. W„ 1. 

Substituting the selected first approximation into (19'), we derive better values for o k, and 
. We determine ' -< such that does not contain a constant term. The expansion of 3g has the 
form 


t CO 

7"' (A^. . A_^ r^'). 


The expansions of fg and are of similar type: 


+ X 



a 


r-2) 




j * 1 


Substituting (35) and (36) into (19'), we obtain the system of "linear" equations: 

4 j2 . 4 (1 . m) j *1 '<o] + I % = A+j' 

’ f -^0 ^ + i 4 [4 j2 _ 4 (1 - m) j - I a.. - A_. . 


(35) 


(36) 


(37) 

(38) 
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We obtain from (37) and (38): 


where 


a+i = tj. 0] A^. + (j. 0) A_. (39) 

a - j = (j, 0) A^. + {j, 0} A_., (40) 


(2j m^) 

[j, 0] = + f 

4 (4 j2 _ , 2 m2) 


^ 4 j 2 (4 j 2 _ + 2 m 2) 

( 2 j *l+m) 2 +i(,<^+ m 2 ) 

tj, 0} . :: [ - j, 0]. 

4 j 2 (4 j 2 _ ^ 2 m2) 

The coefficients [ j, 01, (j, 0) and { j, 0 } remain the same from iteration to iteration. We substitute 
the values of J <, and Vg thus obtained into (29') and derive without any integration. The 
value of the constant is obtained at each iteration from the condition 

"'o r-,- 

We proceed to higher approximations in the same manner, till we reach the convergence. 


THE TERMS DEPENDING UPDN THE FIRST POWER OF THE ECCENTRICITY 

The computation of the first-order effects in e is intimately connected with the computation 
of Cg . We have already derived the system of differential equations 


1 + m)2 +-(•<. 


3 

2^° 


- M 


(28) 


(D2 - < 4. 2 m2) n ^ (w2 _ 1) (29) 

^ 3 m3 (1 . m D-I) [ d i 'fo^ 

- 3 m2 (1 - m Dgd [ (1 ^ Vg) vj . 
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m2 r^l. 




- « w5 (1 




2 


which govern the perturbations proportional to e and are in a form suitable for solving by iteration. 
The relation 


2 Ye + Cl -■ Ce ^ Cl + ig) 0 

can be used to determine without any additional integration. This equation can serve for checking 
the computations and for deriving the value of . At each iteration step, the differential equations 
for terms with the first power of e as a factor have the form: 



1 







(41) 


CD2 


2 m 2 ) V ^ 




(42) 


and correspondingly, the solutions will be of a similar form 




r 


V, = V (e:. 


Ye = 2^ (^+2i + c 


e f"2‘- 

'=+i 




(43) 


Substituting these expressions into the differential equations and comparing the terms in both sides, 
we derive 


e. = [i, c] + (i, c) /3J., (44) 

= (i, c) + {i, c} /S^., 
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where we set 


(2 i -i- Cg - 1 - m)2 X ^ m2) 

[i. C] = 1 _ ( 45 ) 

(2 i Cp)2 1(2 i + Cg)^ - + 2 m2] 


(i, c) = - ^ ? , (46) 

^ (2 i * c^)2 [(2 i ^ Cg)2 - + 2 m2] 

(2 i i Cjj + 1 4- m)2 * i ( /<^ m2) 

li. cj = f _• (47) 

(2 i 4- Cg)2 [(2 i + Cg)2 _ * 2 m2] 


At each cycle, the value of is obtained from the condition 


c2 r K - 2 m2 



The term o-^ / ^ is of the third order in m. The largest portion of it is contributed by the evection 
in ^ and rj . The remaining, much smaller, part is contributed by smaller terms in a and „ . 
Their order in m is not smaller than four. The deriving of and requires, as in the classical 
approach, a special device, because a spurious small divisor appears in (45) - (47) for i = 0. 

In the method based on the iterative procedure, the values of and can be derived from 
the equations 



1 4- m)2 4 - (k 



(48) 


|(Cg _ 1 - m)2 + ^ +m2)J 





(49) 


which follow directly from (41), and combined with the additional condition 


"o - "d = 1- 


(50) 


which we can borrow from the Hill- Brown exposition. In fact, it is equivalent to the definition 
of the mean orbital eccentricity of the satellite. Equations (48) - (50) are compatible and their 
determinant equals zero. We deduce from the last three equations 


0 = (^0 +1 '"o)/ + 1 - + 2 Kp +im2j, 
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*=0 = (^0 - 1 - "1)2 + 2 /<o + i m2]. 

The substitution of these values into (50) can serve as a check of the speed of convergence of the 
iterative process, and of the accuracy of computations. In applying the present theory to the motion 
of the seventh satellite of Jupiter, we found that the iterative procedure is rapidly convergent, 
even when we start with such a crude approximation as 


" \ 4 4 32 32 J 

f i r' - T - — m C2-" - — m r . 
' \ 4 4 32 32 2 


Only a few iterations are necessary to obtain.- , r. , v c„ accurate to 10'® , which seems to 

be the standard accuracy in the expansion of the general theories of satellites of outer planets. 

If the correct value of is initially given, then coefficients (45) - (47) remain the same from 
iteration to iteration. We can dispense with the differential equation for , and the iterative 
process is considerably shortened. The correct value of can be obtained using Andoyer's (1926) 
equation, which is equivalent to Hill's equation but provides a fast convergence for and is con- 
venient for programming. We follow Hill in the computation of auxiliary quantities, and set: 


. n . D) I ^ (1 - .fg)), 


■■ . (1 - 

D) 

' r-' 

' (1 - V- 


A' -■ (1 - a)"' D 

'1 » 

3' 

- A' (1 - 

•)-> D 

P’ - 1 

1 

O 

(A' 

- B') -r m, 



Q' - i (A' - B'). 

S - - Ck w3 . m2) - (2 P'2- 0' 2 - D O') 





2j 


^2i. 
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Introducing the notations 


*^0 = ^0 * "o' 


Uj = (j ^ + 2 j Cg r rig)-‘, 

we write the Andoyer expansion in the form 


"o " ^2 <^^+2 * ^-2> 


- " f^-2 ^-4'> 

+ 2 (u^2 ^+4 * ^-2 ^~4 ' ^ + 2 ^-2^ 


S2 


S ® (/■‘+2 ^ +4 ^^+6 " ^-2 '^-4 '^-e * ^*2 ^'+4 * ''-2 '^- 4 ^ 


2 ( /, ^2 "+4 '^+6 ' *"-2 '"-4 '^-6 * '■ +2 "+4 ^^+6 


-2 ^-4 ^-6 " '■‘+2 “+4 ■' '■'-2 '■'-4 


' +2 '^-2 +4 ' '^+2 '■‘-2 “- 4 ^ 


+ 2 "^+6 


-‘+2 '^+4 ^+6 ^ ^-2 

'■'-4 '■‘-e * ^ + 2 '^-2 

f ^+4 

^4 '^+6 " ''-6 

2 ,_2 * 2 

M _2 

*' - 2^-6 * '^+ 2'^ +4 

' M _2 ^-4 ^ + 2 “ 

-2 ^ 

2 

'■‘-4 ‘ " +2 “-2 ‘ " 

■ 42 --V 


, ^ " +4 ‘■‘+6 * '' -4 '' 

-6 * *" + 2 '' +6 * -2 

^-0 


■-J 


+ 2 "^-4 " ''-2 '"'+4^ * ^6 ^'‘+6 ‘'-6^ 


The Andoyer expansion is convenient for programming and for computing and by iteration. 
It can be easily extended (Andoyer, 1926), if necessary. We tested it on Jupiter VII and found that 
four iterations were necessary to derive the value of accurate to 10 ® . 


HIGHER ORDER TERMS IN i AND^i 


, -j + • and Zj + ^, 


At this point we change the notations and replace c , ■ and z by .f j + 
respectively. As before, ^ and designate the totality of effects of the first order. 


'■'1 ■ ''e ■ '-c' * 

\ ' ”e' " 
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and ^ , V and s from now on designate the effects of higher orders, starting from the second, 

^ = ^2 + ^3 + ■ ■ ■ ’ V = V 2 * V 3 > S = Z - Zj = Zj + Z3 + . . . . 

From (27) we obtain 

f Dg . 1 ^ m) 2 ^ i + m2)] - I = M ij . N r,, - ( 

where - {^fi/'dri)^ are the first-order terms in - They contain either e' or a as a factor. Taking 

into account 

r = -0, = <f. 

=*6 *e 

we have for the effects of higher orders 




1 - m)^ 


. 3 

2 0 


77 r M ^ N 


( 51 ) 


- 2 J c (Djj + 1 * m) A 


where 


F = - + _ [2 (D . 1 . m) b D ^ (S D) 2 ] ^ (S c )2 + S. ( 52 ) 

d 77 \d 77/^ ^ 

This last expression represents the cross-action of terms of lower characteristics as compared 
to the characteristics of terms to be determined. We now introduce averaging operator w, 
which separates the constant terms from an expansion. Thus 


9K 




0. 0. 0 . O' 


Multiplying both sides of (51) by rj^ and both sides of (28') by g and subtracting the results, we have 

Dq ( Dg ^ 5 Dg 7,,) + 2 (1 + m) (7), D ^ ^ Dg T,J (53) 


(I '<0 - n) ^ F - 2 7,^ S c (D, + 1 . m) A 


Taking into account (43) we deduce easily that 

a«{Do (V, % ? - 7,^)} n 0, 
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n^) =0, 


(54) 


Applying the averaging operator to both sides of (53) and using (54), we deduce the basic equation 
for the determination of - c ; 


gK( F) - 2 5 c g» ; r,^ (Dg * 1 ^ m-) A f J -- 0. (55) 

The computation of terms . in r c is closely associated with the computation of 

terms with the characteristics of the type e,\^, where is a monomial in basic parameters. These 
characteristics produce the term c ' ^ of order , ^ in the expansion of c . Designating by 


F 



■' -r 2 i * C 


a’ . 

~ J 


r- 2i "C 


1 


the aggregate of terms with the characteristic e and with the exponents +(2i + c), we obtain 
from (55) 


C - _ ' ’ -- . 

- 2 . 1 . AfJ 

The terms whose characteristics are different from the type mentioned above do not contribute to 
the expansion of -c. Only terms of g and r, with characteristics of the form e\}, where ' is a 
monomial in the basic parameters, contribute terms of the order of > ^ to the expansion of c. We 
obtain the differential equations for the terms of higher orders and higher characteristics in the 
developments of t and v by expanding F, as given by Equation (52), in powers of the basic param- 
eters. Substituting 

S c = . . 

S D . - D, - . . . 

f - ^3 + • • • 

V = Vj * V3 + • ■ • 

y - Xi ^ X2 ^ ■ ■ ■ 

and using (26), (33) - (34'), we can obtain the expansion of 

- [2 (Dg + 1 + m) S D + (S D)2] ^ (5c)^ +S 


22 



in F without difficulty. The development of the first term 


Bn /B n\ 

B 77 \B 77/^ 


requires the application of the operator 


- exp 



+ 77 



+ z 



in its expanded form to - where ^,-q and z are replaced by , rj^ and z^, respectively. 

The value of is set equal to zero after the differentiation. 

The expansion of r has the same form as T in (31) with 




k 



S 


- 1, 2. 3, . . . , 

and Equations (32) and (32') can be used to obtain the expansions of the corresponding t . .-operators 
in terms of 


From the aggregate of terms of a given order we must extract the terms with a given characteristic. 
Such an extraction operator is imbedded into analytical programming languages. In FORMAC, 
for example, it is associated with the COEFF-symbol. A typical differential equation which we 
have to solve has the form 


(Dg . 1 - . ni^)] ^ --0 ■ M,f * N77 + 


(56) 


where cd are the cross-effect terms with the characteristic • . We may start the iteration by 

X 

neglecting and N 77 and repeat the iterative process till we reach the numerical convergence. 
The typical differential equation we have to deal with at each iteration step has the form 


[(D„ - 1 - . i - m2)] C - I -0 r, - A (A 


A’ , 
— 1 


r- 2i 
S 


'-b 


Its solution is 



'■), 


(56') 


(56") 


where 


r li, a\ * (i> o') Alj. 


(57) 
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(58) 


I. j = ( i , a) Aj + { i , a} AJ j , 


and where we set 


(2 i 4 a - 1 - m)2 4 — (k^ 4 

[i, a] = ? (59) 

(2 i 4 o.)^ [(2 i t - Kjj * 2 m^] 


(i, a) = -^ (60) 

2 (2 i + a)2 [(2 i + a)2 - t 2 m^] 


(2 i - a 4 1 ^ h (k^ t m^) 

{i,a} (61) 

(2 i .. a)^ [(2 i 4 a)2 _ Kj, 4 2 m2] 


The imposition of condition (50) has an interesting and important implication. If the character- 
istic has the forme ^ then the right side of (56) and the corresponding parts of and contain 
the terms and . As a consequence of (50) the coefficients ;y and will be equal (Brown, 
1896 and 1899) ° 


If we use (57) - (61) to determine and , we shall have small numbers in the numerators and 
denominators, and consequent loss of accuracy. Knowledge of (62) resolves the difficulty associated 
with the appearance of a spurious small divisor. From (56') and (56") we obtain at each iteration 
step 


[(Cg 4 1 - 4 i * m2)] 4 I 4:^ 4 A^, 

[(Co-l-m)2 .m2)] v„=A;. 


and, taking (62) into account, we deduce 


^^0 ^ -^0 - Vt ^ 1 - m)2 . 2 K, 4 i m2], (63) 

^0 " ^0 " ^’o - 1 - m)2 4 2 Kp 4 i m2] ; (64) 
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should be determined from (63). Equation (64), because of the relative smallness of the de- 
nominator, gives a less accurate result, but can be used for a check. 

PERTURBATIONS IN THE Z-COORDINATE 

We transform the differential equation (3) to the form suitable for computation of z and by 
the process of iteration. We have from (3) 

(■Dg-m^-/<w3')z = - [2D|j SDi-(SD)2_^(w3-vv3)] ■ (65) 


This form of the differential equation for z emphasizes the role of the operator 2 Dq SD in formation 
of the linear equation for terms of the expansion of g. The term 

- l2D|j SD + (bD)^-/<(w^-Wg)] z 

contains cross- effects. Using the substitution 


Z = T Ul + - (I ^ ^6] + 

which resembles the substitution in the theory by Andoyer (1926), and taking (21) and (22) into 
account, we obtain instead of (65) 


(D^ - - k) - K (w* - 1') \p 

+ 2 7 (1 + m - gg) D .fo + r® D Vq) 

+ r (/<o + m2 - g2) [(1 + - (1 + v^) T"] 

+ 1 m2 7 [(1 + 77g) - (1 + - 2 Dg <> D Zj + Z, 

where we set 


z = - (S D)2 Zj 


[2 Dg S D + (S D-)2] (z 


Zj)-k(w2_Wq^)z 


1 B n 

2 B z 


(67’) 


or, taking into account (17) and 


P2 


Z 


1 ’ 


Z = - Zj (S g)2 _ [2 Dg S D + (S D)2](z - z^) _ k (3 w= y + B) 


B n 

”3 z 


(67") 
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Equation (67) is simpler than the corresponding equation of Andoyer, It has the form which empha- 
sizes the relative importance of different terms in \jj and, because of the presence of the small 
parameter m^ in - 1, it is suitable for solving by the iterative process. Setting 

0-= Zj ^3 + • • ■ 

where is of the first order in y and cr is the aggregate of all terms of higher order and higher 
characteristics. We can split equation (67) into two equations, 

(D^ _ - «) i/'j - K (w^ _ 1) (68) 

+ 2 7 (1 + m - . r® Dg 7,3 ) 

+ 7CK3 + m2 - g2) [(^1 + ^3) J+e _ (1 + ^3) 

+ |m2 > [(1 , 7,3) - (1 - e'o^ 

and 

(D2 - m2 - k) ct = « (w^ _ 1) s - 2 - 2 5g D3 r Zj. (69) 

The coefficient of the typical term 

^+g + 2j _ ^"28'2j_ j y_ J 

in the expansion of the right side of the equation (68) will be of the order of m 2 1 ' ! . Thus these terms 
will decrease very rapidly. The term 


^+g-2 _ ^-g+2 

constitutes an exception, because it will be affected by a small divisor in the process of integration. 
Originally of the order of m^, this term will be of the order of m after the integration is performed, 
and will be the largest term in the expansion of 4 ,^ . The value of g3 is obtained from the condition 
that no term of the form 


should be present . Substituting (35) and (36) into (68), and designating by Kg the coefficient 
at the term 


_ j-g 

in the expansion of k (w^ - 1) 0^ , we have by suppressing this term in the right side of (68); 
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= ^<0 +m2 +|m2 a_j + Jfg. 

In the last relation, only K,, varies from iteration to iteration. This method of computing can 
replace a method based on the use of the determinantal equation or its analytic equivalent. 
Another method we suggest is the application of Andoyer's expansion to Adams' equation (1877): 


+ 03 



— aj 


Substituting the expansions of .f and rj into (68) and keeping only the most influential term, we 
obtain an approximate differential equation for : 


J-D2 _ m2 - «• (l -| a2, -| a_2jj 

>{[(l - 2m X |m2) -(go -2)2] a.^ + | m^} ( ^ - 2 _ . x ; 


The solution of this differential equation 


[(l - 2 m . I m2) - (g^ - 2)2] a.j 


3 2 
2 


(go - 2)2 - m2 - -f “ l^-i) 


(/;+e-2 _ r,'K + 2-) 


(70) 


serves to start the iterative process. At each iteration step, equation (68) has the form 


(D2 _ m2 - ,<) i/Zj 3 ^]] C^+b+ 2, _ ^-g-2j-) 


j X 1, 2, 3 M. --- 0 (m2j), 


and the solution is 




k. (^ + g + 2j _ 


(71) 
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where 


k. = M. / [(gg + 2 _ k] . (72) 

If gp is known in advance, then the only term in the right-hand side of (68) which changes from 
iteration to iteration is 


" _ 1 -) 4,^. 

The other terms, and the denominator in (72), remain the same from iteration to iteration. 
Multiplying equation (69) by and equation 

(D2 ^ - fr V/3-) = 0 

by s , and subtracting the results, we obtain 

°o ‘^^1 *^0 z^ Z - 2 i g Dg r Zj. (73) 

Taking into account 

5DJ<D^ (z, D(,a - aD^, z,)} 0, 

we deduce from Equation (73) 


®Uzj Z) - 2 ' g TO(z, '■ z,) ■ 0- 


(74) 


or 


■ (1 ■ - (1 - '/o') Z> 

- 2 g 5K(z, r Zj) 0. 


(75) 


Equation (74) or (75) serves for the determination of og. Evidently only these parts of F having 
characteristics of the form > • ^, where is a monomial in e , e' , y and - 1 . , contribute to the 
expansion of g. They produce terms of the order of ' z in hg . 

We derive the differential equations for the terms of higher orders and higher characteristics 
in cr by expanding z in powers of the basic parameters. 


Substituting 


g ■ 

- gj 

* ^4 + 

D 

D, 

- - 


2 

4 

V 

Vl 

\*2 
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Z = Zj +Z2+Zj + '''‘ 

into (67") and taking (34) and (34') into account, we obtain the development of the first three terms 
in (67") without any great difficulty. The development of the fourth term is derived by applying 
the operator T, in its expanded form, to - 1/2 3n/Bz in which f , -n and z are replaced by , 

Vq andZg = 0, respectively. 
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CONCLUSION 


We discussed certain features which, it seems to us, are peculiar to a form of the theory of 
satellites in rectangular coordinates, if an analytic programming language is used. 

The works of Hill and Brown can hardly be surpassed for mathematical elegance and for deep 
insight into the form of the expansion of coordinates and frequencies and their dependence upon the 
constants of integration. 

Here Hill and Brown have said the last word, providing trigonometrical series are used as the 
means of solution. Space is left, however, for the decision as to howto derive the solution and what 
the most convenient form of the differential equations would be if the analytic manipulations are 
performed automatically. 

In our approach to these problems, we found inspiration in the work of Andoyer on the Hill-Brown 
lunar theory. We found that it contains several interesting points, which sound quite modern. This 
beautiful work deserves revived attention. 

We obtained Hill's periodic intermediate orbit by solving simultaneously the differential 
equation for the complex coordinate and the Andoyer differential equation for the parallax, bypassing 
the formation of the infinite system equations for coefficients of the variation terms. 

We transformed the differential equations of the Hill-Brown theory to a form which contains a 
small parameter in the right-hand sides and thus favors the application of the method of iteration 
and the use of analytical programming language. 

Together with the expansion of the coordinates, we derive the expansion of the parallax. Its 
knowledge facilitates the expansion of the right-hand sides of the differential equations in powers of 
the constants of integration. The right hand sides of the differential equations contain T. . the 
polynomials in , X2 • • • • The representation of T. ^ in the form of the polynomials in the con- 
stants of integration constitutes an essential part of the theory and programming. The knowledge 
of this representation facilitates greatly the formation of the differential equations for the terms 
with the prescribed characteristic. 

We expanded the differential operators into series in squares of the constants of integration. 

Only the terms of order zero in the expansions of the differential operators are employed in 
the integration of the differential equations. The remaining terms are responsible for producing 
the cross effects between the perturbations of different orders. 

All programming was done by R. Baxter, to whom the author is greatly indebted for his 
generous help and friendly cooperation. He is using FORMAC as a basis. 

We are aware that FORMAC is not the most efficient language as far as computer time is con- 
cerned. However, it has been valuable to us in terms of flexibility and ease of programming of the 
numerous sub-problems into which lunar theory in coordinates is split. It was also very useful in 
checking the correctness and the applicability of the theoretical thought. 
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We are applying the theory presented here to the seventh satellite of Jupiter and compare the 

results with the work of Tokmalaieva (1956). 

Goddard Space Flight Center 

National Aeronautics and Space Administration 
Greenbelt, Maryland, September 28, 1970 
188-41-51-06-51 
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